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A STANDARD KINEMATIC 
MODEL FOR FLIGHT SIM- 
ULATION AT NASA-AMES 

by 

Richard E. McFarland* 


SUMMARY 


The use of a standard kinematic model provides for 
effective utilization of the flight simulation facilities at 
NASA-Ames. This paper describes the model as well as the 
functional relationships used in its derivation. 


♦Programs Manager, Computer Sciences Corporation, Mountain View, 
California. 
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1 . 0 Introducti on 


The material in this paper constitutes the kinematic model 
used in the generalized aircraft simulation program structure, 
herein called BASIC, and additionally outlines certain stan- 
dardized auxiliary simulation components such as the atmospheric 
and turbulence models. The kinematic model is common to all 
of the simulation computers at NASA-Ames, and is maintained as 
is a general computer library. Because of the diverse problem 
areas this model is called upon to investigate, it is extensive. 

In the following development, the equations presented 
should not be construed as constituting either the computa- 
tional load or the procedural order of the digital programs 
comprising BASIC. The actual computations, which accommodate 
multi-loop, real-time programmi ng structures, are segmented 
into modules which are executed at rates based upon their 
relative frequency content. 

The various coordinate frames used in the model have the 
advantage of great relative accuracy without the necessity of 
resorting to computationally expensive techniques, such as 
double precision arithmetic. Although integrations are per- 
formed in a pseudo-inertial frame, with compass-directional 
axes, the mathematical relationships are rigorous, and resul- 
tant positional information may be used for certain naviga- 
tional studies. 

All axes systems used are orthogonal, right-handed triads. 
2.0 Earth (E) and Local (L) Frame Relationships 

The Earth Frame is the principal frame for the development 
of the mathematical relationships used in BASIC since, neglect- 
ing extraterrestrial considerations, it neither rotates nor 
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translates. However, it is not an especially interesting 
frame in the problem of aircraft simulation, and quantities 
related to this frame do not appear explicitly in any of the 
finalized equations. 

The Earth Frame is an inertial frame with the origin of 
coordinates at the earth’s center (spherical); the axis 
intersects the North Pole, and the axis intersects the 
zero degree longitude line (Greenwich) at "zero time". 

The Local (L) Frame is situated on the earth's surface, 
directly under the vehicle. Its X L axis points northward and 
its axis points eastward; both are parallel to the earth's 
surface. The Z L axis points towards the earth's center. This 
is not an inertial frame since it follows the motion of the 
aircraft. Its distance from the center of the earth (R g ), 
however, is constant. 

The E and L Frames are related through the Inertial 
Longitude Angle (xj) and the Latitude Angle (x) as shown in 
Figure 1. Since the earth rotates about the Z^ axis with the 
angular velocity u , a vehicle's longitude x will differ from 
the Inertial Longitude Angle xj as a function of elapsed time: 

X = X I - U) g t (2.1) 

The inertial latitude does not differ from that of a rotat- 
ing earth. 
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Figure 1 Description of E-L Axes Systems 


The transformation from the Earth Frame to the Local Frame 
(see Appendix A for the derivation) is given by 



- S x C t i 

- S \ S T J 

CA 

^ T E 2 L ^ 

- Sx I 

CT I 

0 


- C x C t j 

- C x S t i 

-S\ 

the trigonometric functions 

"si ne" 

and 


( 2 . 2 ) 


abbreviated "S" and "C" throughout. 


Since the L-Frame is always under the aircraft on the earth's 
surface, the center of the L-Frame is always a distance R g 
(20,898,908 feet) from the center of the E-Frame, which is lo- 
cated at the earth's center. The Z L always points downward 
(perpendicular to a spherical earth's surface). 


Whenever the aircraft translates over the earth's surface 
a relative rotation between the L and E Frames occurs which may 
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be represented by the instantaneous rates (p L , q^, r L ) of the 
L-Frame itself (X L> , Z L ). The omega-cross operator [« L ], 


[n L ] 



(2.3) 


and two identities which relate a transformation matrix to its 
derivative in terms of this operator. 


^E2L^ _ "t n L^ T E2L-* 

(2.4) 

[t L 2 E ] = [ T L2E^lJ 


are developed in Appendix B. The second derivative follows as 


where 


and 


^E2lJ ^ T E2L^ 


[n L ] 


0 - f L q L 
r L 0 -p L 

-q L p l o 


[n L ] 2 


-r 2 -q 2 P L q L r,_p L 

P L q L “ r L~ P L r L q L 

P|_ r L q L r L _q L _p L 


(2.5) 


( 2 . 6 ) 


(2.7) 


- 4 - 



From Appendix A, the instantaneous rates about the L-Frame 
may be related to the translational velocity of the aircraft 
over the earth's surface by 


( 2 . 8 ) 


This vector introduces the Coriolis effect into the model; its 
influence, being velocity proportional, is understandably infin- 
itesimal for most applications. However, problems of even coarse 
navigation aboard an SST would suffer from its absence. 



The inertial position vector (E-Frame) is related to the 
position vector in the Local Frame by 



[ t L2E ] 


Z L“ R e 


( 2 . 9 ) 


where is the negative of the altitude h . Taking into account 
the fact that R g is a constant, the inertial velocity vector is 
the time derivative of the above equation: 



( 2 . 10 ) 
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However, in the L-Frame, no motion occurs from the origin of 
coordinates except along the axis (X^=Y^=0). Defining 

R = R - Z , = R p + h (2.11) 

e L e 


the inertial velocity vector (transformed to the instantaneous 
Local Frame) is 




— 


— — 


— ‘ 

V N 


*E 


0 


0 

V E 

T E2L 

*E 

= 

0 

+ [ & ^ ] 

0 

V p. 


Zr 


-R 


-R 



E 


— — 




which, using (2.3) and (2.8) is also 


. 


*— 


• 

V N 


-Rq L 


RA 

V E 

= 

Rp L 

= 

Rt jCos A 

V D 


-R 


-R 


— 


— 


(2.13) 


From these elements, (2.8) may be rewritten 


p L 

1 

R 

V E 

“ V N 

_ r L_ 


-V E Tan A 


(2.14) 


The inverse radius factor reiterates the argument that this 
vector generally has a very small magnitude. 


The inertial acceleration vector is the time derivative of 
the inertial velocity vector (2.10), 
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x 'e 

Y E 

[f L2E ] 

X L 

y l 

+ 2 r f L2E J 

X L 

y l 

* t T L 2E ] 

X L 

' y l 

Z E 


Z L- R e 


Z L 


Jl_ 


[T 


L2E 


(2.15) 


X L 

y l 

+ 2[n L ] 

X L 

y l 

+{( C« L ] 2 + [0 L ]} 

1 

_1 _J 

X >- 

1 

Z L 


_ Z L_ 


Z. -R 
L e 


and again. 


since motion occurs only along the axis, 


— 

/ 

— 


““ — 


— — i 

X E 

\ 

0 


0 


0 

Y E 

^ T L2E^j 

0 

+ 2[n L ] 

0 

+ { [Q l ] 2 + [fl L ]} 

0 

_ Z E_ 

V 

-R 


R 


_-R_ 


(2.16) 


This is the inertial acceleration vector, and as such, may be 
transformed to the Local Frame and equated to the total acceler- 
ation vector (both applied and field forces considered) 



(2.17) 


where the force of gravity is related to the sea level weight 
W or gravitational acceleration g Q by 

F g = W ( R e /R) 2 = mg 0 (R e /R) 2 (2.18) 
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3 . 0 Local (L) and Body (B) Frame Relationships 


The Body Frame uses conventional aircraft notation: the X g 
axis passes through the nose of the vehicle, the Y g axis points 
toward the right wing, and the Zg axis passes through the bottom 
of the vehicle. The Body Frame origin is located at the vehicle 
center of gravity. 

The transformation matrix from the L-Frame to the B-Frame 
(Tl 2 b) is an identity matrix whenever the aircraft is parallel 
to the earth's surface (horizontal) and pointing northward. 

From Appendix C, 



CeCij/ 

c es^ 

-Se 


II 

1 — 1 
CO 
C\J 
_J 
t— 

t 1 


S<f>SeS^+C 4 >Ci}j 

S<j>Ce 

(3.1) 


C<f>SeC^ + S<|>S^ 

c<f>seSij7-s<j>c^ 

C 4 .Ce 



Forces applied in the B-Frame may be represented in the 
L-Frame by the inverse (or transpose) relationship 


f n 

^ T L2B^ 

—1 

X 

_l 

f e 

f ty 

f d 


_ F tz_ 


and with the addition of the field force this is proportional 
to the inertial acceleration vector (2.16), transformed to the 
L-Frame (see (2.14) and (2.17)) 



_f n “ 


-2q L R-r L P L R-q L R 

1 

m 

f e 

= 

2p L R-r L q L R+p L R 


f d +f g 


-R+(q^+p2)R 


_ 






The derivative of components of (2.14) with substitutions 
from (2.13), 


Pl 

1 

*E + VL 

Jl_ 

R 

3 + V D q L _ 


(3.4) 


when substituted into (3.3), produces the Local Frame force 
vector. This illustrates the individual contributions due to 
coriolis and centripetal accelerations as well as the applied 
and field forces. 



f n 


\ y D- r L?L R+ h 

1 

m 

LU 

U_ 

= 

" P L V D" r L q L R+ ^E 


f d +f g 


( q L + P L ) R + ^D 


(3.5) 


By substitutions from (2.14), this equation is used to obtain 
the time differential of the inertial velocity vector trans- 
formed to the instantaneous Local Frame: 






" q L V D +r L P L R 

^E 

-IE 

II 

f e 

+ 

P L V r L q L R 



_ f d +f g_ 


-(qL + pf )R 


(3.6) 



f n 


V N V D’ V E TanX 

1 

m 

LU 

Li— 

♦J 

V E + n ^ 


f d +f g 


_3 +v e> _ 
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These accelerations are integrated to form the Local Frame 
velocity vector, and this integration includes the effects of 
both the earth's rotation and rotations produced by vehicle 
transl ati on . 


V 



V 


E 


V 


D 



dt 


( 3 . 7 ) 


The second-order Adams-Bashford algorithm is used. The 
vehicle velocity over the earth's surface is as above, minus 
the vector contribution of the earth's rotation (u> ) about the 
axis. The Local Frame velocity of the vehicle with respect 
to the earth's surfa c e is: 


V N 


V N 


r° 

S i n A 0 


0 

V EE 

= 

V E 

“ 0) 

e 

-Si n A 

0 -Cosa 


0 

V D 


_V 


0 

Cosa 0 


R 



' V N~ 



0 




= 

V E 

<D 

3 

I 

Cos A 

i 





V D 



0 




Only the eastward component is influenced. 


Not to be confused with a “north wind", which emanates 
from the north, in BASIC terminology a "velocity north" heads 
northward (if positive), etc. The total wind vector is treated 
as a random variable with both a trend and random component 
represented by 





— — — — 




V TWN 


V NW 


V NTURB 

V TWE 

= 

V EW 

+ 

V ETURB 

TWD 


_yDW_j 


_!!dturb_ 
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in which the random portion may originate in the Body Frame, 
as in 





— 

V NTURB 


U TURB 

V ETURB 

^ T B2L^ 

V TURB 

V DTURB 


W TURB 


( 3 . 10 ) 


The aerodynamic velocity vector is then the difference 
between the vehicle velocity with respect to the earth's surface 
and the total velocity of the air mass with respect to the earth's 
surface 


V NR 


V N 


V TWN 

V ER 

= 

V EE 

- 

V TWE 

V DR 


V D 


V TWD 


In the Body Frame this is expressed as 


U B 


V NR 

V B 

^ T L2B^ 

V ER 

W B 


V DR 


( 3 . 12 ) 


with resultant magnitude 


RW 


= l /u B 2 + V B 2 tU B 2 ' 


=T/vj!+v7F+v 


71 


NR ER 'DR 


( 3 . 13 ) 


which is not to be confused with the ea rth-rel ati ve velocity 
magni tude 


V = Tv' 2 + v 2 + v 21 

V T V N V EE V D 


( 3 . 14 ) 
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or the ground speed 


V 


G 



V *"+V 
v N v 



(3.15) 


The body Axis acceleration vector is the derivative of (3.12), 
(see Appendix B and (3.25)) 


— — 


— — 1 


— 


— 

°B 


0 -r T q T 


U B 


^n'^twn 

^B 

= - 

r T 0 " P T 


V B 

+ t T L2B^ 

V*TWE 



-q T p t o 


W B 


^D'^TWD^ 


(3.16) 


In Section 9.0, however, it is shown that our turbulence model 
develops disturbances in velocity space only, and the differen- 
tiation of random variables to obtain pseudo-acceleration terms 
for (V TWN , V TWE , V TWD ) is not an attractive discrete operation. 
Furthermore, within the BASIC system the vector ( U g , Vg, Wg) is 
an open-loop calculation with the following two exceptions: 

(1) During the i ni tial -condi ti on trim process (see (11.0)) in 
which random components are necessarily eliminated, and (2) in 
the computation of a and B (3.30 and 3.31), which are seldom- 
used variables; the introduction of turbulence into these terms 
should best be left to the discretion of researchers. For these 
reasons the random elements of (3.16) are suppressed by use of 
the approximation 


— 


— — 

'V^TWN 


'V^'nw 

V*twe 

A/ 

y*EW 

V V TWD 


V -V 

V D DW 


(3.17) 


which includes only acceleration terms of the "trend" of the 
wind profile. Finite differences are adequate for this computa- 
ti on. 


-1 2 - 



Equation (3.16) is verified by [1]; notice that the integra- 
tion of (3.16) in the reference frame (B), which may have a high 
rotational velocity, is judiciously avoided by (3.7). This 
approach is basically that developed by L. E. Fogarty and R. M. 
Howe in early publications on "Analog Computer Solution of the 
Orbital Flight Equations". 

The rate of change in longitude over a rotating earth is 
equal to the inertial rate of change due to translation about 
the axis, minus the earth's rate itself (derivative of 2.1) 
such that (2.13) may be rewritten 




— . 



V N 


RX 


0 

V E 

= 

Rt Cos X 

+ Ru> Cosx 
e 

1 

V D 


-R 


0 


— — 


— — 


And by comparison with (3.8), 


H • 

1 

1 

V EE /Cosx 

— 
1 

1 

*3| 

_ V N 


(3.18) 


(3.19) 


where "t" is the vehicle longitude and "a" is the latitude. A 
direct integration, with appropriate initial conditions, produces 
the position vector 


T 



X 




(3.20) 


which determines the position of the vehicle with respect to 
the surface of a spherical, rotating earth. This integration 
is performed with the modified Euler algorithm. 
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By assuming that the X-Z plane is a plane of symmetry, the 
Euler moment equations [2] may be written 

L T = PB I XX +q B r B^ I ZZ" I YY^'^ r B +P B q B^ *XZ 

M T = q B I YY +P B r B^ I XX _I ZZ^'^ r B _P B^ *XZ (3.21) 

N T = r B I ZZ +P B q B^ I YY' I XX^"^ P B" q B r B^ ! XZ 


By simultaneous solution the angular acceleration couple in 
these equations may be eliminated. The resultant form of Euler's 
Dynamical Equations is 



( c 1 r B+ C 2 p B ) q B 

5 r B P B + C 6 ^ r B -p B ^ 
( C 8 P B +C 9 r B ^ q B 



in which the coefficients have been computed in the following 
order: 
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C 0 



C 1 


C 0^ I YY" I ZZ^ I ZZ' I Xz} 

C 2 


( iO I XZ^ I XX -I YY + I ZZ^ 

C 3 


c o I zz 

C 4 


c o ! xz 

C 7 


1 Y Y 

C 5 


C 7^ZZ _I XX^ 

C 6 


C 7 1 XZ 

C 8 


C o{^ I XX _I YY^ I xx +I xz} 

C 9 


C 0 I XZ^ I YY" I ZZ" I XX^ 

1 

o 

o 

1 


C I 
l 0 XX 


The moments of inertia, and hence these coefficients, are 
usually designated as constant for simulations in which there 
is no concern with real-time variations in either mass or 
center- of- gravity. This fact, however , does not restrict BASIC 
to either the selected principal axes or the stationary coef- 
ficient assumptions. BASIC's highly modularized form permits 
easy access to (and replacement of) groups of equations such 
as these. 
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The Body-Frame angular rates are gained through the inte- 
gration of (3.21) by use of the second-order Adams-Bashford 
algorithm. 


PB 

°*B 



(3.24) 


The rotation of the L-Frame, using (2.14), is referenced to 
the B-Frame through the transformation 






P LB 


p L 

q LB 

= [ T L2B^ 

q L 

r LB 


_ r L_ 


(3.25) 


and the difference between this rate vector and that resulting 
from the applied torques in the B-Frame produces the differences 
in the orientation of the Local and Body Frames. 


> 


~p B 


Plb 

q T 

= 

q B 

- 

q LB 

r T 


r B_ 


r LB 


(3.26) 


This equation is therefore used to produce the angular rate 
differences which in turn are used to create the Euler sequence 
rates (see Appendix C). The Euler angle sequence must be speci- 
fied - in BASIC it is assumed (in transforming from the Local 
Frame to the Body Frame) that the sequence order is (1) yaw, 

(2) pitch, and then (3) roll. 
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$ 


(qySinf+ryCos<}))/Cose 

• 

e 

= 

q^Cos^-rySi n<j> 

♦ 


p^+ij Sine 


The Euler angles which relate the B-Frame to the L-Frame, 
in the specified sequence, are computed from the above vector 
by use of the modified Euler algorithm. 



( 3 . 28 ) 


Special programming considerations permit the simulation 
of hovering and rearward flight. The angles of attack and 
sideslip are 


Tan _1 [W B /U B ] (-|-<a<^) 

(3.29) 

Tan' 1 [V B /(j”l/u B 2 + W B 2 *) ] (-ir<B<ir) 

(3.30) 


in which j is unity, with the sign of Ug. The rates of change 
in these quantities are 


U B W B" W B U B 

u b 2 *“b 2 


(3.31 ) 


e = 


[(U b 2 +W b 2 )V b -V b (U b U b+ W b W b )] (3.32) 


' 2 1/u 2 tU 2 ' 
RW » U B W B 
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4.0 Velocity Vector Initialization 


The total velocity with respect to the earth's surface may 
be described in terms of pointing angles which have no aero- 
dynamic relationship whatsoever: 


V N 


CoSYyCoSY^j 

V EE 

= V T 

CosYySi ny H 

_ V D _ 


-SinY v 


Velocity initialization in BASIC is the simultaneous satis- 
faction of the constraints of (4.1) regarding the direction of 
the total velocity vector, and the selected aerodynamic velocity 
magnitude given by one of the two options: 

( V -a (option #1 ) 

V RW = T (4- 2 ) 

^ V Q /( . 59248 5 "1/p/ p 0 ') (option #2) 

With the first option, V Q is used as an input for Mach 
Number which is multiplied by the speed of sound (at altitude). 
With the second option V Q is used for the initial Equivalent 
Airspeed in the units of knots; combining the appropriate con- 
version factor and the indicated function of atmospheric density 
ratio, it too results in the aerodynamic velocity V RW 

The equivalent airspeed is 

V eq . .592485 V R „y^7V < 4 ' 3 ) 

It should be noticed that the aerodynamic velocity magnitude 
is that which is specified regardless of the magnitude of direc- 
tion of an arbitrary wind profile. In addition, the initializa- 
tion procedure also requires the specification of the direction 
(y h and y v ) of the total velocity vector. 


- 18 - 



Initialization of a vehicle on the ground is handled 
differently. Equation (4.1) is solved directly with the 
assumption that the taxi speed is the input V Q (it replaces 


5 . 0 Pilot and Center of Gravity Relationships 

The position and acceleration of both the pilot's station 
and the center of gravity are computed in BASIC for reasons 
which include the control of visual and motion simulators. 

The Body-Frame inertial acceleration of the center of 
gravity (which would be sensed by accelerometers placed there) 
relates to the applied forces proportionally 


Q) ^ 

X 

1 

1 

m 

f tx 

3 Y 

f ty 

a Z 


f tz 


such that at the pilot's station the sensed accelerations are 



(5.2) 
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where, as in (2.3) et al 


[« B ] 


0 " r B q B 
r B 0 " P B 

“ q B P B 0 


(5.3) 


The position vector (X p , Y p , Z p ) relates the pilot station 
to the center of gravity in the Body Frame; its derivatives are 
general ly zero . 


By specifying the runway parameters longitude "ip", latitude 
"Xp", height above sea level "hp" and clockwise direction with 
respect to North "e R ", the radius from the center of the earth 
to the runway is 


r R ■ R e + h R (5 ' 4 > 

and the center of gravity of the vehicle is located north, east 
and above the runway according to 


ANr 


X-Xp 


0 


aE r 

- r r 

(x-Tp)CosXp 

+ 

0 

(5.5) 

h CG 


0 


_ h “ h R_ 



In terms of the Runway Frame, which is defined with x-axis 
"down" the runway and y-axis "to the right", with origin at the 
threshold, the center of gravity position is further amplified: 


X CG 


Cose R 

S i n e p 


ANr 

y cg 


- S i n e p 

Cos 0 R 


aE r 

__ 


— 

— 





(5.6) 
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The pilot's station is related to the runway by 


aN pr 


J 

PD 

1 

+ ^ T L2 B ^ 

X P 

aE pr 

= 

aE r 

Y P 

_h pr 


" h CG 


_ Z P_ 


X PR 


Cose R 

S i n e R 


aN pr 

y pr 


- S i n 9 R 

Cose R 


^ e pr_ 


Equations (5.5) through (5.8) are "flat-earth" approximations 
which are only valid in the vicinity of the runway. 

6. 0 Atmospheric Quantities 

BASIC provides many typical atmos phe r i c- rel ated quantities; 
some examples appear below. 

The total temperature ratio is related to the mach number 

T R = 1 + .2M 2 (6.1 ) 

If the velocity is subsonic, the total pressure ratio is 

P R = (T r ) 7/2 (6.2) 

but if it is supersonic it is 

P R = (166.9 )M 2 / (7-M~ 2 ) 5/2 (6.3) 


These two expressions are equal if M=1 . Below h = 36,089 
feet the ambient temperature and pressure ratios are given by 


T ar = 1 -6 . 875X1 0" 6 h 
P AR = T AR 5 - 256 


(6.4) 


- 21 - 



Above this altitude they are 


T ar = .751895 

P = . 2 234e [ " 4 - 806X10 " 5(h ' 36089)] 


(6.5) 


The ambient temperature (°C) has provision for an incremental 
temperature, and is given by 

T ft = aT a + T Ap (518.69)/ (1 .8) (6.6) 


The ambient pressure is 


P fl - (2116.2 )P ar 


The total temperature and pressure are 


T 

P 


T 

T 


T R T A 


P R P A 


and the impact pressure is 

q C = P T' P A 


(6.7) 


( 6 . 8 ) 


(6.9) 


BASIC utilizes the 1962 ARDC tables for atmospheric density 
(p) and the speed of sound (a), both as functions of altitude 
(h) up to 240,000 feet, with data points every 2000 feet. 
Optionally, a constant relationship may be selected. 

The atmospheric density and speed of sound may be modified 
by del ta-temperature effects as follows 
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p 


p TABLE< T A- AT A>/ T A 


( 6 . 10 ) 


a " a TABLE _ V T A / ^ T A" AT A^ ' 

The mach number and dynamic pressure are 

M = V RW /a 

q = Y p v rw 

and the calibrated air speed is computed from 

V c = (.592485)a 0 “l/ 5 |- (1+ - c/2116 2)- 2857 -l] * 


( 6 . 11 ) 


( 6 . 12 ) 


where "a 0 " is the speed of sound at sea level. 
7 . 0 Forces and Moments 


The atmospheric and other quantities shown above are 
generally useful to BASIC users for utilization by their Engine 
and Aerodynamic modules. These modules, together with the total 
landing gear influence, are then used to produce the total 
applied forces and moments used in (3.2) and (3.21) respectively. 

The total force, applied at the vehicle center of gravity 
and with all components transformed to the Body Frame, consists 
of aerodynamic, engine, and gear contributions (discussed in 
Section 8.0): 


f tx 


f ax 


f ex 


f gx 

f ty 

= 

f ay 

+ 

f ey 

+ 

f gy 

f tz 


f az_ 


_ f ez_ 


f gz_ 
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Contributions to the moments about the vehicle center of 
gravity in the Body Frame are computed from the same sources 
as are the forces: 


l t 


l a 


l e 


l g 

m t 

= 

m a 

+ 

m e 

+ 

m g 

n t 




n e_ 


_v 


(7.2) 


This over-simplified structure is designed only for modular 
isolation of components, especially in program-debug mode; addi- 
tional contributions to forces and moments are of course per- 
mitted. 

8 . 0 Landing Gear 

The user is provided (from BASIC) the values of compression 

and stroke rate of each gear and is expected to return to BASIC 

values for oleo force (F. ), friction (F D ), and side force (F<~ ). 

L i K i 5 i 

The modeling of these forces is the user's responsibility. They 

usually involve considerations of static loading, damping, braking 

and steering. 

The height of a selected position near the tail of the air- 
craft is also monitored 

h^. = h^g + X-j-Sine - Z^Cose (8.1) 

to determine whether or not the ground has been intercepted. 

If so, the vertical velocity of intersection is 

hy = fi + e (XyCose + ZjSine) (8.2) 
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Each landing gear is monitored by use of various elements 
in the [T^g] arra y of Appendix C 


h G i ~ ^ h CG“ T 13 X G i “ T 23 Y G i " T 33 Z G i ^ /T 33 


(8.3) 


A negative component denotes that the particular gear is 
on the ground. Its rate of compression is computed with a 


difference equation. 

Also, 

a program "flag" is set 


( 0 

Gear not on ground 

\ 

■ 

(8.4) 

( 1 

Gear on ground 


for event marker indication, landing gear lights, etc. 


With the above information the user's Gear Routine provides 
the necessary forces for the equations 


F RX i ~ F R i “ F L i 6 

F RY . = F L + F S . (8.5) 

F D7 = F, + F d Sine - F c Sin<j> 

RZ i L i R i S i 


Addi ti onal ly 
(F T p) may be 


a tail reaction force (Fy^) and friction force 
developed from (8.1) and (8.2). 


The gear contributions to the total moments, including 
tail impact, are 
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L G =E 


F RZ i Y G i ‘ F RY i (h G i + Z G 


A 


’g = " X T F TR + E[ F RX i (h G. + Z G i ) ' F RZ i X G i 


( 8 . 6 ) 


N G = E 


F RY i X G i “ F RX 




The total force contributions due to the landing gear and 
tail are 


— — 


■ 


— — 

f gx 


f tf 


F RX i 

f gy 

= 

0 


F RY i 

f gz 


f tr 


_ FRZ i_ 


The subscript "i" above runs from one to three which corre- 
sponds to the nose, right and left wheels respectively. 


9.0 T u r bulence 


As in the case of aerodynamic velocity turbulence (3.10), 
computations are also made for angular velocity disturbances. 

The resultant variables are not, however, used in the angular 
acceleration computations, but rather, are simply made available 
to the user for inclusion in his aerodynamic coefficient buildup 
equations [3]. 


P BWN 


P B 


P TURB 

q BWN 

= 

q B 

+ 

q TURB 

r BWN 


r B 


r TURB 


(9.1 ) 
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Figure (2) is a Laplace diagram of how four separate, digital 
white noise sources (Z^) may be filtered to produce aircraft dis- 
turbances conforming to the Dryden spectra model. The sample 
data correction is shown along with the continuous equivalent of 
what is actually simulated using Z-Transforms with a zero-order 
hold. Parameters are supplied by the user, and they consist of 
the wing span (b), the characteristic lengths (Ly, Ly, L w ), and 
the dispersions (a^, ay, o^). 

Our digital model is based upon the original work of Messrs. 
Neuman and Foster [4] of NASA, and has since gained applicability 
to a wide range of aircraft configurations. Certain mainframe 
manufacturers have their own versions of this model, notably 
Boeing; their analysis concludes that they have produced filters 
which "match the Von Karman PSD 1 s for m < 1 0 rad/sec". This model 
is also available with BASIC. It is noteworthy that the tech- 
niques utilized call for the introduction of turbulence effects 
in the Earth Frame rather than the Body Frame. 


10.0 Integration 


The equations of state given in (3.7), (3.20), (3.23) and 
(3.27), are integrated in BASIC by use of the second-order Adams- 
Bashford predictor and the modified Euler algorithm which are 
given by 


y n+ l = y n + (3y n -y n _ 1 ) At/2 ( Adams-Bashford 2nd) 
y n+l = y n + ( y n +y n-l ) At / 2 (modified Euler) 


( 10 . 1 ) 


The selected integration algorithms (and their interaction) 
have undergone an extensive phase and magnitude analysis (in the 
total system environment) over forcing function frequencies of 
up to a few hertz. This was accomplished by the comparison of 
non-real -time segments computed with unity-millisecond cycle 
times to segments obtained by use of realistic cycle times in 
the tens-of-mi 1 1 i seconds range. 
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Arbitrary transfer functions are usually simulated with 
constant-cycl e-time state space techniques such as Z-Transf orms . 

A general technique, which we call SOLACE, has the advantage that 
on-line changes may be made to the coefficient structure and 
system form. This system, which utilizes matrix exponentiation 
techniques with either a zero or first-order hold, will generally 
give the real-time solution to the following stationary system [5] 


y 


B m + l S ” 

+ B S m-1 + • • • 
m 

+ B l" 

C n + l S " 

, r c n_ l , 

t c n s * ••• 

+ C 1 


( 10 . 2 ) 


with a computation overhead equivalent to the multiplication of 
an mxn-matrix and an n-vector. 


11 . 0 Trimming Capability 

A routine which drives the state vector (p B , q B> r B , U B , V B , 
Wg) to zero by appropriate manipulation of an arbitrary set of 
controls (up to six) is fully integrated into the BASIC structure. 
The routine uses the classical techniques of regression analysis 
and acts as a forcing function to the entire simulation rather 
than to subsets thereof. This process is indicated in Figure 3. 
Automatic features of the system include a process which senses 
that the vehicle is on the ground, ignores the normal control 
list, and performs a weight-on-wheels trim which takes the gear 
geometry into account. Asymmetrical flight conditions do not 
pose any problems to our 6-deg rees-of- freedom trim routine 
(called BQUIET) providing that the steady state solution exists 
within the aircraft's data map and control deflection extrema. 

12.0 Auxiliary Relationships 

The standard aircraft simulation model includes many "canned" 
relationships which accrue to each facility user. They result 
in instrument drives, indicator lights, event markers, strip 
chart displays, localizer and glide slope signals, and other 
computations which include 

(a) Outer and middle marker events resulting from the vehicle 
being within a vertical cone emanating from selected 
pos i ti ons , 
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FIGURE 3 - TRIM-SIMULATION INTERFACE 
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(b) Altitude "trips", resultant from the wheels being 
below 

(1 ) 1500 feet 

(2) 200 feet 

(3) The flare height (pre-set constant) 

(4) A selected height (pilot variable) 

(c) Thrust asymmetry computations resulting in event 
marker indicators, 

(d) Automatic strip chart speed control as a function of 
altitude and computer mode, and automatic strip chart 
run number identification and calibration after each 
run, and 

(e) Landing gear transit delays and discretes. 

13.0 Concluding Remarks 

In order to utilize BASIC it is generally required that 
the user furnish the information necessary to program modules 
(subprogams) which "fill in the blanks". This minimally con- 
sists of functional curves and tabulations, geometrical consid- 
erations, schematics for the control system and engine charac- 
teristics, aerodynamic coefficient buildup equations, other 
data, and at least the gear static loading curves. Then, in 
a pseudo-closed-shop environment, programmer-analysts ascertain 
the correct frequency domain division of labor, allocate special 
COMMON symbols, program the modules, link them to all BASIC 
COMMON symbols, interface the program to the simulation facility, 
perform static and dynamic checks, and assist the project engi- 
neers in the gathering of data, its interpretation and reduction, 
and in the expansion and elaboration of the simulation capability. 
It is, of course, an integral part of this process that the 
researcn goals be published at project initiation. 
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A non-real -time version of BASIC is available on NASA-Ames' 
IBM-360 computer for the purpose of off-line debug procedures 
and all-digital data acquisition and analysis. By the standard- 
ized use of "plug-in’ 1 modules containing the specific aircraft 
configuration information, programs developed on this computer 
are directly applicable (in source form) to the real-time simu- 
lation computers. Figure 4 outlines the batch-processi ng-type 
executive routine which is used to supplant the real-time execu- 
tive controller in these applications. 

Appendix D is included to acquaint prospective users of the 
NASA-Ames real-time flight simulation facility with our data 
handling and analysis capability. 

Complete COMMON lists for the BASIC system are available upon 
request. Over 600 variables are involved, sorted both numerically 
and alphabetically. 
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figure: 4 

BATCH PROCESSING EXECUTIVE 
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Appendix A 

E-L Frame Relationships 


The transformation from the Earth Frame to the Local Frame 
involves two rotations, in the following order: 

[ t e2L ] = (1) Yaw " T i"’ (2) pitch T 1 ' X " 



Where "roll" is defined as a positive rotation about the 
x-axis, "pitch" is defined as a positive rotation about the 
y-axis, and yaw is defined as a positive rotation about the 

z-axi s . 
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Tne derivative of ( abov e) is 

- C x C t j x + S a S t j t j - C a S x jA-SxCt jT j -Sax 
CijTj - S x j x j 0 

S A Ct j A + C A S t j t i S A S x i A - C A C x j x j -CAA 

and from Equation (2.4) it is 


r i S t i + q lC A C x i r^Cx j +q^CASx j q^SA 

r l S A C x j - p ^ C A C x j r S A S x j - p^C A Sx j -r^CA-p^SA 

q^SxCxj + p^STj - q^ S A S x j - p^Cx j q ^ C A 


From an examination of the elements. 
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Appendix B 


Matrix Operations 


The vector cross product wxv has an equivalent form in 
matrix notation which may be given by [ft]v. To show this [6] 
we take two vectors with the same basis (i, j, k) 


w = iw + jw w + kw 

a y z 

V = iv x + jv y + kv z 


and perform the indicated operations 

i j k 


u = wxv = 


WWW 

x y z 

v v v 
x y z 


i (w yV z -w 2 Vy) + J(w z v x -w x v z ) 
+ ^ W x V y- W yV x ) 


w V -W V 


V 

y Z Z y 


X 

w zV w x v z 

■ [fl ] 

w J 

v y 

> 

3 : 

> 

3 


v 

l_x y y x 1 


z 


(B-l) 


(B-2) 


By inspection of each term, the omega-cross operator (ex- 
pressed as elements of the leading cross-product vector) must 
be the skew-symmetric matrix 


0 


-w 


z 


w 


y 


[%]= 


W 


2 


0 -w x 



(B-3) 
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Or, if elements of the second vector are used. 


5 - [o v ] 


w 


X 



(B-4) 


where again by inspection. 


[“,] ■ 




( B - 5 ) 


Note the change of sign when elements of the trailing vector are 
used. 


The transpose of (B-3), for instance, shows that a post- 
operator may also be used: 


u 


[u v O = ([ Q ]v) = * C fl ] 

s y l 


Tv V V 1 

L x y z J 


-V [a] = - 


•w_ w. 


0 


-w. 


-"y "x 0 


(B-6) 


Consider some vector v in an arbitrary orthonormal coordinate 
system undergoing simultaneous rotations about all three axes. 

Let p be the instantaneous rotational velocity about the x-axis; 
q and r are similar quantities defined with respect to the y and 
z axes. Define the vector w = (p, q, r). 


Consider first a rotation about the x-axis (roll) through 
an angle <f>. The vector coordinates of some point P change to Q 
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due to the rotation 


x = x 

y ' = y cos 4 > + z si n«j> (B-7) 

I 

z = z cos<j> - y sin4> 

or if <J> = pAt, with At an infinitesimal time, 

x ' = x 

y = y + z p A t (B-8) 

z = z - ypAt 

Likewise, a rotation about the y-axis (pitch) in time At gives 

I 

x = x - zqAt 

y' - y < B - 9 > 

» 

z = z + xqAt 

and a rotation about the z-axis (yaw) in time At gives 

I 

x = x + yrAt 

y=y-xrAt (B-10) 

I 

z = z 

If we simultaneously perform all three rotations, and note 
that the non-commutati ve terms are second order in At [7], we 
see that 

x' = x + (yr-zq ) At 

y = y + (zp-xr)At ( B - 1 1 ) 

z' = z + (xq-yp)At 
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At 


X -X 

y -y 

i 

Z -Z 


X 

y 

z 




(B-l 2 ) 


Hence, the time rate of change of a vector undergoing rotation 
only (not change in length) may be described by the cross product 
of (B-3), where w represents the instantaneous angular rates of 
the axis undergoing rotation. By simple vector partitioning [8], 


“11 

“12 

“13 

“21 

“22 

“23 

“31 

“32 

“33 



"11 

“1 2 

“13 

“21 

“22 

“23 

“31 

“32 

“33 


(B-l 3 ) 


whenever the length of each column is preserved. This is true 
for any orthogonal matrix, e.g. the transformation from point P 
to point Q, so that 


T PQ = “ I-^q] T pQ 


(B-l 4) 
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Appendix C 

L-B Frame Relationships 


The rotational sequence in transforming from the Local 
Frame to the Body Frame is 


[T L2B ] = (1) Yaw (2) Pitch "e", (3) Roll *V 


1 0 0 

0 C<(> S</> 

0 - S <j> C (j) 

CeCi/i 

C 4> S 9 Ci// + S</i Si// 


Ce 0 -Se Ci// 
0 1 0 -Si// 

Se 0 Ce 0 

CeSij/ -Se 

S</)SeSi/j + C<t>Ci// S C e 

C (j/S eSi//- S<j)Ci// C 4> C 6 


Si// 

Ci// 

0 


(c-1) 


From (2.4) and (3.26) this transformation is related to the 
rates which produce its variations by the relationship 

[t L2B ] = - [a T ] [t L2B ] (c- 2 ) 

which may be expanded (dispensing with superfluous subscripts) 



*11 

*1 2 

*13 


0 

r T 

■‘»T 


T 11 

T 1 2 

T 1 3 



*21 

T 22 

T 23 

= 

- r T 

0 

p T 


T 21 

T 22 

T 23 

( C-3 ) 


*31 

T 32 

T 33 


_ q T 

-Pt 

0 


-J 31 

T 32 

T 33_ ; 


Selecting 

only 

three 

of these 

ni ne 

relationships 

and 

equati ng 


them to the derivative of (C-1) results in 
T -j -j = - 0 S 0 C^“ tjfC 0 Sl|) 
r T T 2 1 “ q T T 31 

= r-|.(S<f>SeCij/-C<f>Si/i)-q-|.(C<f)SeCi/i + S<)>Si//) 


( C - 4 ) 
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- 0 C 0 


r T T 23“ q T T 33 
= r-j.S<j>Ce-qjC<i)Ce 


T 


33 


-<t>S<j>C 6 + 0 C 4> S 9 

q T T l 3 ” p T T 23 
- q jSe- p-pS 4 > C 6 


(C-5) 


( C-6 ) 


These three equations, when solved simultaneously, yield 
equation (3.27). The selection of a different Euler rotational 
sequence would produce a different form of this angular rate 
equation. 
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Appendix D 

Data Handling and Analysis 


Normally, 24 digi tal -to-anal og converters are reserved for 
on-line, strip chart recordings of 48 multiplexed data channels; 
the resultant time histories, plus postmortem printouts and plots, 
often constitute the entire data acquisition phase - and this is 
usually sufficient for analyses such as work-load determination 
or handling quality studies. However, for any sophisticated 
analysis of simulation data, our real-time, digital, magnetic 
tape, data acquisition system should be used; we call this system 
"RUNDUM". 

RUNDUM makes it possible to record up to 

N = 15. At (limited by 1000 due to buffer size) 

floating point variables (where "At" is the frame time, in milli- 
seconds, at which output occurs). The resultant magnetic tape(s) 
produced are in a format which is compatible with all of the 
laboratory's data handling and analysis programs. 

The data handling programs available are varied, flexible 
and optimized for data-transf er rate, i.e. the computational 
stream is limited by the mechanics of the tape handlers. For 
instance, this permits the " RUNDUM- to- I BM 360 conversion" pro- 
gram to transfer data at the rate of about 8000, 32-bit words 
per second. 

Currently available data handling routines include: 

(1) "WRITE360" - Converts a RUNDUM tape to an IBM-360- 
compatible tape. 
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(2) "STRIP" - Makes strip chart recordings from selected 
RUNDUM segments. 

(3) "XPLOT" - Cross-plots selected RUNDUM data. 

(4) "MACPRINT" - Lists selected data on a line printer. 

(5) "TTEST" - Tests and sizes a RUNDUM tape by displaying 
(a) the variable names, (b) the maximum and minimum 
of each variable, and (c) variable values every 200 
points (frames). It performs this operation for each 
file (simulation run) on the tape. 

The laboratory's primary tool for data analysis is MAC/RAN* 
[9], a system which has been created for the analysis of time- 
series data. This system operates directly on data from a 
RUNDUM tape. 

The MAC/RAN system provides the following capabilities to 
simulation experimenters: 

(1) Filtering - A wide range of filter types and algorithms 
are available prior to data analysis. 

(2) Statistics - Means, variances, higher order moments, 
tests for randomness, histograms, and other operations. 

(3) Auto-correlation and cross-correlation functions. 

(4) Power spectral density and cross spectral density 
functions. 

(5) Coherence functions. 


♦MAC/RAN is a registered trademark of Measurement Analysis 
Corporation, Marina del Rey, California. 
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